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A variational theory is developed to study electrolyte solutions, composed of interacting point-like 
ions in a solvent, in the presence of dielectric discontinuities and charges at the boundaries. Three 
important and non-linear electrostatic effects induced by these interfaces are taken into account: 
surface charge induced electrostatic field, solvation energies due to the ionic cloud, and image charge 
repulsion. Our variational equations thus go beyond the mean-field theory, or weak coupling limit, 
where thermal fluctuations overcome electrostatic correlations, and allows one to reach the opposite 
strong coupling limit, where electrostatic interactions induced by interfaces dominate. The influence 
of salt concentration, ion valency, dielectric jumps, and surface charge is studied in two geometries, 
i) A single neutral dielectric interface (e.g. air-water or electrolyte-membrane) with an asymmetric 
electrolyte. A charge separation and thus an electrostatic field gets established due to the different 
image charge repulsions for coions and counterions. Both charge distributions and surface tension are 
computed and compared to previous approximate calculations. For symmetric electrolyte solutions 
close to a charged surface, two zones are characterized. In the first one, in contact with the surface 
and with size proportional to the logarithm of the coupling parameter, strong image forces and 
strong coupling impose a total ion exclusion, while in the second zone the mean-field approach 
applies, ii) A symmetric electrolyte confined between two dielectric interfaces as a simple model 
of ion rejection from nanopores in membranes. The competition between image charge repulsion 
and attraction of counterions by the membrane charge is studied. For small surface charge, the 
counterion partition coefficient decreases with increasing pore size up to a critical pore size, contrary 
to neutral membranes. For larger pore sizes, the whole system behaves like a neutral pore. For 
strong coupling and small pore size, coion exclusion is total and the counterion partition coefficient 
is solely determined by global electroneutrality. A quantitative comparison is made with a previous 
approach, where image and surface charge effects were smeared out in the pore. It is shown that 
the variational method allows one to go beyond the constant Donnan potential approximation, with 
deviations stronger at high ion concentrations or small pore sizes. The prediction of the variational 
method is also compared with MC simulations and a good agreement is observed. 

PACS numbers: 03.50.De,87.16.D-,68.15.+c 



I. INTRODUCTION 

The first experimental evidence for the enhancement 
of the surface tension of inorganic salt solutions com- 
pared to that of pure water was obtained more than eight 
decades ago [TJ |2] ■ Wagner proposed the correct physi- 
cal picture [3] by relating this effect to image forces that 
originate from the dielectric discontinuity and act on ions 
close to the water-air interface. He also correctly pointed 
out the fundamental importance of the ionic screening of 
image forces and formulated a theoretical description of 
the problem by establishing a differential equation for the 
electrostatic potential and solving it numerically to com- 
pute the surface tension. Using series expansions, On- 
sager and Samaras found the celebrated limiting law ^ 
that relates the surface tension of symmetric electrolytes 
to the bulk electrolyte density at low salt concentration. 
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However, it is known that the consideration of charge 
asymmetry leads to a technical complication. Indeed, im- 
age charge repulsion, whose amplitude is proportional to 
the square of ion valency, leads to a split of concentration 
profiles for ions of different charge, which in turn causes 
a local violation of the electroneutrality and induces an 
electrostatic field close to a neutral dielectric interface. 
Bravina derived five decades ago a Poisson-Boltzmann 
type of equation for this field [5] and used several ap- 
proximations in order to derive integral expressions for 
the charge distribution and the surface tension. 

These image charge forces play also a key role in slit- 
like nanopores which are model systems for studying ion 
rejection and nanofiltration by porous membranes (see 
the review [3] and references therein, and [7] for a re- 
view of nano-nuidics) . Several results have been found in 
this geometry and also for cylindrical nanopores beyond 
the mean-field approach (using the Debye closure and the 
BBGKY hierarchical equations) and averaging all dielec- 
tric and charge effects over the pore cross section. Within 
these two approximations, the salt reflection coefficient 
has been studied as a function of the pore size, the bulk 
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salt concentration and the pore surface charge. 

More precisely, the strength of electrostatic correla- 
tions of ions in the presence of charged interfaces without 
dielectric discontinuity is quantified by one unique cou- 
pling parameter 

3 = 2^i 2 B a s (1) 

where q is the ion valency, and a s the fixed surface 
charge [STUOj. The Bjerrum length in water for mono- 
valent ions, £b — e 2 / (^T:e w kBT) s» 0.7 nm [e w is the 
dielectric permittivity of water) is defined as the dis- 
tance at which the electrostatic interaction between two 
elementary charges is equal to the thermal energy k B T. 
The second characteristic length is the Gouy-Chapman 
length £q = l/(2nq£ B a s ) defined as the distance at which 
the electrostatic interaction between a single ion and a 
charged interface is equal to k B T '. The coupling param- 
eter can be reexpressed in terms of these two lengths 
as 3 = q 2 £ B /£a- On the one hand, the limit 5^0, 
called the weak coupling (WC) limit, is where the physics 
of the Coulomb system is governed by the mean-field 
or Poisson-Boltzmann (PB) theory, and thermal fluctu- 
ations overcome electrostatic interactions. It describes 
systems characterized by a high temperature, low ion va- 
lency or weak surface charge. On the other hand, 5 — > oo 
is the strong coupling (SC) limit, corresponding to low 
temperature, high valency of mobile ions or strong sur- 
face charge. In this limit, ion-charged surface interac- 
tions control the ion distribution perpendicularly to the 
interface. For single interface and slab geometries, several 
perturbative approaches going beyond the WC limit [TT] 
or below the SC limit US] have been developed. 
Although these calculations were able to capture impor- 
tant phenomena such as charge renormalization |14j , ion 
specific effects at the water-air interface |T51 [TB] , Man- 
ning condensation [17] , effect of monopoles [18] or attrac- 
tion between similarly charged objects, they also showed 
slow convergence properties, which indicates the inabil- 
ity of high-order expansions to explore the intermediate 
regime, 5 ~ 1. This is quite frustrating since the com- 
mon experimental situation usually corresponds to the 
range 0.1 < 5 < 10 where neither WC nor SC theory is 
totally adequate. 

Consequently, a non-perturbative approach valid for 
the whole range of 3 is needed. A first important attempt 
in this direction has been made by Netz and Orland [T5] 
who derived variational equations within the primitive 
model for point-like ions and solved them at the mean- 
field level in order to illustrate charge renormalization ef- 
fect. Interestingly, these differential equations are equiv- 
alent to the closure equations established in the context 
of electrolytes in nanopores [6]. They are too compli- 
cated to be solved analytically or even numerically for 
general 3. A few years later, Curtis and Lue [5D] and 
Hatlo et al. [21j investigated the partition of symmetric 
electrolytes at neutral dielectric surfaces using a similar 
variational approach (see also the review [H]). They have 
also recently proposed a new variational scheme based on 



a hybrid low fugacity and mean- field expansion |23j , and 
showed that their approach agrees well with Monte-Carlo 
simulation results for the counterions-only case. How- 
ever, this method is quite difficult to handle, and one has 
to solve two coupled variational equations, i.e. a sixth 
order differential equation for the external potential to- 
gether with a second algebraic equation. Within this 
approach, these authors generalized the study of ion-ion 
correlations for counterions close to a charged dielectric 
interface, first done by Netz in the WC and SC limits [24], 
to intermediate values of 3. They also studied an elec- 
trolyte between two charged surfaces without dielectric 
discontinuities at the pore boundary, in two cases: coun- 
terions only and added salt, handled at the mean-field 
level [25] . Although this simplification allows one to focus 
exclusively on ion-ion correlations induced by the surface 
charge, the dielectric discontinuity can not be discarded 
in synthetic or biological membranes. Indeed, it is known 
that image forces play a crucial in ion filtration mecha- 
nisms [B]. The main goal of this work is to propose a 
variational analysis which is simple enough to intuitively 
illustrate ionic exclusion in slit pores, by focusing on the 
competition between image charge repulsion and surface 
charge interaction. Moreover, our approach allows us to 
connect nanofiltration studies [26-28 with field-theoretic 
approaches of confined electrolyte solutions within a gen- 
eralized Onsager-Samaras approximation 4J character- 
ized by a uniform variational screening length. This vari- 
ational parameter takes into account the interaction with 
both image charge and surface charge. We also compare 
the prediction of the variational theory with Monte-Carlo 
simulations [25] and show that the agreement is good. 

The paper is organized as follows. The variational 
formalism for Coulombic systems in the presence of di- 
electric discontinuities is introduced in Section ^] Sec- 
tion |III| deals with a single interface. We show that the 
introduction of simple variational potentials allows one to 
fully account for the physics of asymmetric electrolytes 
at dielectric interfaces (e.g. water-air, liquid-liquid and 
liquid-solid interfaces, see ref. [30]), first studied by Brav- 
ina [5J using several approximations, as well as the case of 
charged surfaces. In Section [TV] the variational approach 
is applied to a symmetric electrolyte confined between 
two dielectric surfaces in order to investigate the prob- 
lem of ion rejection from membrane nanopores. Using 
restricted variational potentials, we show that due to the 
interplay between image charge repulsion and direct elec- 
trostatic interaction with the charged surface, the ionic 
partition coefficient has a non-monotonic behaviour as a 
function of pore size. 



II. VARIATIONAL CALCULATION 

In this section, the field theoretic variational approach 
for many body systems composed of point-like ions in the 
presence of dielectric interfaces is presented. Since the 
field theoretic formalism as well as the first order varia- 
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tional scheme have already been introduced in previous 
works p~9| . 120 ) . we only illustrate the general lines. 

The grand-canonical partition function of p ion species 
in a liquid of spatially varying dielectric constant e(r) is 



^ n w.iv 

i=l Nl =0 iVl!A * 



3Ni 



3=1 



,-(H-E s 



(2) 



where At is the thermal wavelength of an ion, pi denotes 
the chemical potential and Ni the total number of ions of 
type i. For sake of simplicity, all energies are expressed 
in units of ksT. The electrostatic interaction is 



H=- J dr'drp c (r)ti c (r,r')p c (r / ) 



where p c is the charge distribution (in units of e) 



P Ni 
<=1 3=1 



(3) 



(4) 



and Qi denotes the valency of each species, p s (v) stands 
for the fixed charge distribution and v c (r, r') is the 
Coulomb potential whose inverse is defined as 



k R T 



V [e(r)V<5(r - r' 



where e(r) is a spatially varying permittivity. The self- 
energy of mobile ions, which is subtracted from the total 
electrostatic energy, is 



1 = 1 



(6) 



where v h c (r) = is/r is the bare Coulomb potential for 
e(r) = e w . After performing a Hubbard-Stratonovitch 
transformation and the summation over Ni in Eq. ([2]), 
the grand-canonical partition function takes the form of 
a functional integral over an imaginary electrostatic aux- 
iliary field <p(r), Z — J T><p 



The Hamiltonian is 



H[<p] 



dr 



i 

where a rescaled fugacity 



8n£ B {r) 



ip s (r)<f>(r) 



(8) 



has been introduced. The variational method consists in 
optimizing the first order cumulant 



F v — F + (H — Hq)o. 



(9) 



where averages (• • • )o ar e to be evaluated with respect to 
the most general Gaussian Hamiltonian [T5] . 

M<t>\ = \ J [m-^)]%\ry)W)-i<t>o{*')\ 

(10) 



and Fo = — ^trlnuo- The variational principle consists 
in looking for the optimal choices of the electrostatic 
kernel v Q (r,r') and the average electrostatic potential 
4>o(r) which extremize the variational grand potential 
The variational equations <5i 7 ! u /i5v ( j" 1 (r, r') = 
F v /S(j)o(r) — 0, for a symmetric electrolyte and 
„ yield 

A^o(r) - 87r^s?Ae-^ w Wsmh[# (r)] 

= -4*e B p.(r) (11) 

-Av (r,r') + 87r£ s g 2 Ae-^ w «cosh[g^ (r)]i;o(r,r') 

= 47r£ B (5(r-r'). (12) 




where we have defined 



W(r) = lim [u (r, r') - v b c (r - r')] 



(13) 



whose physical signification will be given below. The 



second terms on the LHS of Eq. (Ill and of Eq. (12) 



have simple physical interpretations: the former is Att£b 
times the local ionic charge density and the latter is 
Airtsq 2 times the local ionic concentration. The rela- 



(5) tions Eqs. (fill) — (|12[) are respectively similar in form to 



the non-linear Poisson-Boltzmann (NLPB) and Debye- 
Hiickel (DH) equations, except that the charge and salt 
sources due to mobile ions are replaced by their local 
values according to the Boltzmann distribution. On the 
one hand, Eq. (11) is a Poisson-Boltzmann like equa- 



tion where appears the local charge density proportional 
to sinh^o- This equation handles the asymmetry in- 
duced by the surface through the electrostatic potential 
cj>o, which ensures electroneutrality. This asymmetry may 
be due to the effect of the surfa ce cha rge on anion- and 
cation-distributions (see Section III B ) or due to dielec- 



tric boundaries and image charges at neutral interfaces, 
which give rise to interactions proportional to q 2 , and in- 
duce a local non-zero </>o for asymmetric electrolytes (see 
Section III A). On the other hand, the generalized DH 



(7) equation Eq. (12), where appears the local ionic concen- 



tration proportional to cosh^o, fixes the Green's function 
«o(r, r') evaluated at r with the charge source located at 
r' and takes into account dielectric jumps at boundaries. 

These variational equations were first obtained within 
the variational method by Netz and Orland [19]. They 
were also derived in Ref. [21] within the Debye closure 
approach and the BBGKY hierarchic chain. Yaroshchuk 
obtained an approximate solution of the closure equa- 
tions for confined electrolyte systems in order to study 
ion exclusion from membranes [B] . 



Equations ( 11 )— ( 12 ) enclose the limiting cases of WC 



(S — > 0) and SC (S — > oo). To see that, it is interesting to 
rewrite theses equations by renormalizing all lengths and 
the fixed charge density, /5 s (r), by the Gouy-Chapman 
length according to f = r/£ G , p s (r) = e G p s (r)/a s (a s 
is the average surface charge density). By introducing 
a new electrostatic potential <^o( r ) = 90o( r )7 on e can 
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express the same set of equations in an adimensional form 

A^ (f) - Ae"5 sinh0 o (f) = -2p s {v) (14) 

-Au (r, f) + Ke-^ W(i) cosh 4>o(t)v (t, f) 

= 47T(5(f - f')(15) 



where vo — vo£g/£b, W = W£g/£b and we have also 
introduced the rescaled fugacity A = 8tt\£ g 'E. [32] ■ Now, 
one can check that, in both limits S — > and S — > oo, the 
coupling between </>o and in Eq. (11) disappears and 



the theory becomes integrable. Finally, it is important to 
note that this adimensional form of variational equations 
allows one to focus on the role of i>o(r, r) whose strength 
is controlled by H in Eqs. (14) and (Tl5|). However, even 



at the numerical level, their explicit coupling does not 
allow for exact solutions for general H. 

In the present work, we make a restricted choice 
for vo(r,r') and replace the local salt concentration in 
the form of a local Debye-Hiickel parameter (or inverse 
screening length) /-c(r) in Eq. (12), 



n{vf = 8Tre B q 2 Xe-^ w ^ cosh [#o(r)] , (16) 

by a constant piecewise one n v (r) — k v in the presence 
of ions and K v (r) — in the salt- free parts of the sys- 
tem. Note that it has been recently shown that many 
thermodynamic properties of electrolytes are successfully 
described with a Debye-Hiickel kernel [33] . 

The inverse kernel (or the Green's function) v (r,r') 
is then taken to be the solution to a generalized Debye- 
Hiickel (DH) equation 



-V(e(r)V) + e(r)^(r)] v (r,r') 



knT 



tf(r-r') (17) 



below, its explicit form depends on the confinement ge- 
ometry of the electrolyte system as well as on the form 
of e(r). The variational equation for the electrostatic 
potential [31] 5F v /5(j)o(z) = yields regardless of the 
confinement geometry 



d 2 <f> 



At: I, 



^W{z)- qi 4, a (z) 



The second-order differential equation (21), which is 



simply the generalization of Eq. (11) for a general elec- 



0. 



(21) 



trolyte in a planar geometry, does not have closed-form 
solutions for spatially variable W(z). In what follows, 
we optimize the variational grand potential F v using re- 
stricted forms for the electrostatic potential (f>o(z) and 
compare the result to the numerical solution of Eq. (21 ) 
for single interfaces and slit-like pores. 
The single ion concentration is given by 



Pl (z) = A 4 e-Tr^W-^W 



and its spatial integral by 



dzpi(z) = — A. 



dF v 



(22) 



(23) 



We define the Potential of Mean Force (PMF) of ions of 
type i, $i(z), as 



= - In 



Pb 



By defining 



w{z) = W{z) - W b 



(24) 



(25) 



with the boundary conditions associated with the dielec- 
tric discontinuities of the system 



lim w (r, r') = lim v (r,r'), 



(18) 



lim e(r)Vv (r,r') = lim e(r)Vv (r, r') (19) 

r— >£ r— + 

where S denotes the dielectric interfaces. We now re- 
strict ourselves to planar geometries. We split the grand 
potential ^ into three parts, F v = F\ + F2 + F3 , where 
Fi is the mean electrostatic potential contribution, 



S / dz { - 



p s (z)4> (z) 



■|ff( 2 )-((rf„W 



(20) 



F 2 the kernel part and F3 the unscreened Van der Waals 
contribution. The explicit forms of F2 and F3 are re- 
ported in Appendix [A] The first variational equation is 
given by dF v /dn v = d (Fi + F 2 ) /dn v — 0. This equa- 
tion is the restricted case of Eq. (12). As we will see 



where Wb is the value of W{z) in the bulk and comparing 
Eqs. (p2l-p4l, we find 



-w b 



In 7? = Mi-ln(p 6 A?) 



(26) 
(27) 



Hence, qfWb/2 is nothing else but the excess chemical 
potential of ion i in the bulk and q?W{z)/2 = In 7,(2) 
is its generalization for ion i at distance z from the 
interface. They are related to the activity coefficients 7* 
and Ji(z). Note that the zero of the chemical potential is 
fixed by the condition that (f> vanishes in the bulk. The 



PMF, Eq. (26), is thus the mean free energy per ion (or 



chemical potential) needed to bring an ion from the bulk 
at infinity to the point at distance z from the interface, 
taking into account correlations with the surrounding 
ionic cloud. 

Before applying the variational procedure to single and 
double interfaces, let us consider the variational approach 
in the bulk. In this case, the variational potential 4>o is 
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FIG. 1: (color online) Geometry for a single dielectric inter- 
face (e.g. water-air) (a) and double interfaces or slit-like pores 
(b). 
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equal to 0, and the variational grand potential F v only 
depends on k v . Two minima appear: one metastable 
minimum k° at low values of k v , and a global minimum 
at infinity (F v — > — oo for k v — > oo) which is unphysical 
since at these large concentration values, finite size effects 
should be taken into account. It has been shown by intro- 
ducing a cutoff at small distances [5D], that, for physical 
temperatures, this instability disappears and the global 
minimum of F v is k„ ~ Kb given by the Debye-Hiickel 
limiting law: 



Hi = ln(p 6 Af) - \n h l B 
k\ = Mb J2i QiPi,b 



(28) 



From Eq. (27), we thus find Wb ~ —k^Ib and the poten- 



tial w(z) reduces to 

w(z) w v (z, z) - v b c (Q) + K b £i 



(29) 



which will be adopted in the rest of the paper. Further- 
more, problems due to the formation of ion pairs do not 
enter at the level of the variational approach we have 
adopted. Let us also report the following conversion re- 
lations 



h 

K h 



0.l9(K b £ B ) 2 mol.L" 
3.29\/7 run -1 , 



(30) 



where I = ^/2J2i1iPi,b is the ionic strength expressed 
in mol.L - . Finally, the single-ion densities are given by 



Pi (z) = Pitb e-^ w ^-^ z l (31) 



III. SINGLE INTERFACE 

The single interfacial system considered in this Section 
consists of a planar interface separating a salt-free left 
half-space from a right half-space filled up with an elec- 
trolyte solution of different species (Fig. [I]). In the gen- 
eral case, the dielectric permittivity of the two half spaces 
may be different (we note e the permittivity in the salt- 
free part). The Green's function, which is chosen to be 
solution of the DH equation with e(z) — e9(—z) + e w 0(z) 



FIG. 2: (color online) Potential w(z) in units of fc_gT for e = 
(black solid curve), Eq. (361, and e — e w (red dashed curve) 
and Kb(.B = 4. 



and k(z) = k v 6(z) where 0{z) stands for the Heaviside 
distribution, is given for z > by [S] 



W(z) = £ B (Kb - K v ) 



(32) 



kdk 



where 



A(x) = 



e w yjx 2 + 1 - i 
e w Vx 2 + 1 + i 



and F 2 [Eq. (A4)] can be analytically computed 

.3 



where 



A = A(x -)• oo) = 



(33) 



(34) 



(35) 



The first term on the rhs. of Eq. (34) is simply the volu 



mic Debye free energy associated with a hypothetic bulk 
with a Debye inverse screening length k v and the second 
term on the rhs. involves interfacial effects, including the 
dielectric jump A, and k v . 

For the single interface system, as seen in Section|TTJ F 3 
is independent of k v and fioiz), which means that it does 
not contribute to the variational equations. By minimiz- 
ing Eqs. (20) and (34) with respect to n v for fixed </>o(r) 
and taking V oo, one exactly find the same variational 
equation for for the bulk case. Hence, as discussed 
above, we have k v = Kb given by Eq. (28) and the first 
term of the rhs. of Eq. (32) vanishes. This result was 
obtained in [2D] for the special case <p(r) = 0. It is of 
course not surprising to end up with the same result for 
finite <fi(r) since we know that the electrostatic potential 
should vanish in the bulk. This potential combines in an 
intricate way both the image charge and solvation contri- 
butions due to the presence of the interface. The image 
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force corresponds to the interaction of a given ion with 
the polarized charges at the interface and is equivalent 
to the interaction of the charged ion with its image lo- 
cated at the other side of the dielectric surface. As it 
is well known, the image charge interaction is repulsive 
for e < e w (e.g. water-air interface) and attractive for 
e > e w (the case for an electrolyte-metal interface) [55] , 
The interfacial reduction in solvation arises because an 
ion always prefers to be screened by other ions in order 
to reduce its free energy. Hence, it is attracted towards 
areas where the ion density is maximum (at least at not 
too high concentrations for which steric repulsion may 
predominate). This term is non-zero even for e = e w 
since for an ion close to the interface, there is a "hole" 
of screening ions in the salt-free region (where k v = 0). 
Although our choice of homogeneous variational inverse 
screening length allows us to handle the deformation of 
ionic atmospheres near interfaces that are impermeable 
to ions, it does not allow us to treat in detail the local 
variations in ion solvation free energy arising from ion- 
ion correlations (except in an average way in confined 
geometries where differ from the bulk value of the 

inverse screening length, see Section IV below). 

Equation ( 32 1 simplifies in three cases : 
1) For e — (A = 1), where the solvation effect vanishes 
because the lines of forces are totally excluded from the 
air region [35] , Eq. ( 32 ) reduces to 

-1n b z 



w (z) 



2z 



(36) 



This is the case where the image charge repulsion is the 
strongest (see Fig. [2]) . 

2)A slightly better approximation for e 7^ can be ob- 
tained by artificially allowing salt to be present in the air 
region. This gives rise to the "undistorted ionic atmo- 
sphere" approximation [BJ, for which w(z) in Eq. (36) is 
multiplied by A: 

-2n b z 



W 



z)=A£ B - 



2z 



(37) 



Solvation effects are now absent and salt exclusion arises 
solely from dielectric repulsion. Eq. (371 is exact for ar- 
bitrary Kf, and A = 1, or arbitrary A and Kb = 0. 
3)In the absence of a dielectric discontinuity e = e w (A = 
0), the potential can be expressed as 



w(z) = K b e B f(Kbz) 



(38) 



/(*) = 



[l + xf 



K 2 (2x) 



2x 3 



where K 2 (x) is the Bessel function of the second kind. 
One notices that unlike the case A > 0, the potential has 
a finite value at the interface, i.e. w(0) = Kbis/3- 

We note that in this case of one interface, we have 
liim._j.oo <fio(z) = and the fugacity of each species is 
fixed by its bulk concentration according to 



Pifi = lim pi(z) = Xie~ 



(39) 



A. Neutral dielectric interface 

We investigate in this section the physics of an asym- 
metric electrolyte close to a neutral dielectric interface 
(e.g. water-air, liquid-liquid or liquid-solid interface) 
located at z = (a s = 0). For the sake of simplicity, we 
assume e = 0, which is a very good approximation for 
the air-water interface characterized by e — 1 (see the 
discussion in Ref. [!]). Hence we keep the approxima- 
tion w(z) — wq(z) given by Eq. (36). The electrolyte is 
composed of two species of bulk density p+ and p- and 
charge (q+e),— (<?— e) with q + > g_. In order to satisfy the 
electroneutrality in the bulk, we impose p+q+ = p~q-- 
According to Eq. ( 28 ) , the bulk inverse screening length 



noted Kb is given by 



= ^ni B q_p_ (g_ + 



(40) 



and the variational equation (21 ) for the electrostatic po- 



tential is a modified Poisson-Boltzmann equation 



d 2 c 



dz 2 



+ 4Tre B p ch (z) = 



(41) 



with a local charge concentration 



Pch(z) = p-q- e -^™( z )-1+Mz) _ e —s-w(*)+q-M*) 

(42) 

Equation (41 ) can not be solved analytically. Its numer- 
ical solution, obtained using a 4th order Runge-Kutta 
method, is plotted in Fig. [3|a) for asymmetric elec- 
trolytes with divalent and quadrivalent cations and the 
local charge density is plotted in Fig.J^b). 

Fig. [3] clearly shows that, very close to the dielectric 
interface for z < a, image charge repulsion expulses all 
ions (since p c h(z) ~ exp(— 1/z) has an essential singu- 
larity) and <fio is flat. For z > a, but still close to the 
interface, there is a layer where the electrostatic field is 
almost constant (<f)Q increases linearly), which is created 
by the charge separation of ions of different valency due 
to repulsive image interactions. The intensity of image 
forces increases with the square of ion valency and close 
to the interface, p c h(z) < since we assumed q + > q- 
(the case for Mgl 2 ). To ensure electroneutrality, the lo- 
cal charge then becomes positive when we move away 
from the surface (Fig. |3jb)), and the electrostatic poten- 
tial goes exponentially to zero with a typical relaxation 
constant kj,. Moreover, in Fig. [3^ a) one observes that 
when the charge asymmetry increases, the electrostatic 
potential also increases. Knowing that for symmetric 
electrolytes, <po = 0, our results confirm that the charge 
asymmetry is the source of the electrostatic potential <f>Q. 
Fig. |3^b) is qualitatively similar to Fig. 1 of Bravina who 
had derived an integral solution of Eq. (41 ) by using an 



where we used Eq. ( 22 ) 



approximation valid for Kbts *C 1 [5.J ■ In order to go fur- 
ther in the description of the interfacial distribution of 
ions, we look for a restricted variational function 4>o(z) 
which not only contains a small number of variational 
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(b) 




z/L 




■ electrolyte 1:2 
i electrolyie 1 :4 



z/L 



FIG. 3: (color online) (a) Electrostatic potential <f>o (in 
fcsT u nits) for asymmetric electrolytes: numeric al so lution 
of Eq. (41 1 (symbols) and variational choice, Eq. (44 1 (solid 
lines), for divalent and quadrivalent ions and p- = 0.242 M. 
Variational parameters are ~ 1.4k;,, a/£e = 0.12; 0.21 and 
tp = —0.10; —0.156. (b) Associated local charge density profile 
(thick lines) and anion (dashed lines) and cation (thin solid 
lines) concentrations. 



are constants. Hence, we choose a restricted variational 
piecewise solution 4>o(z) 



00 00 



ip [1 + k^(z - a)] i 



-K<j,(z—a) 



for z < a, 
for z > a. 



(44) 

whose derivation is explained in Appendix [B] The vari- 
ational parameters are the constant potential (p, the de- 
pletion distance a and the inverse screening length re^. 
The grand potential (B5) derived for this solution was 



optimized with respect to the variational parameters us- 
ing the Mathematica software. The restricted variational 
potential (44) is compared to the numerical solution of 
Eq. (41]) in Fig. [3] for electrolytes 1 : 2 and 1 : 4 and 
P-l% — 0.05. The agreement is excellent. One notices 
that the screening of the effective surface charge created 
by dielectric exclusion enters into play when z > Kl ■ Fi- 
nally, let us note that since kjjIb = 1-37 and K\)Ib = 1-77 
respectively for the monovalent and quadrivalent elec- 
trolytes in Fig. [3j the method adopted by Bravina is not 
valid. 

To summarize, the charge separation is taken into ac- 
count by the potential ip (which increases with q + / 'q- ) 
and the relaxation constant ~ 1.4/C& is almost inde- 
pendent of q + /q_. Interestingly, the variational parame- 
ter alls — 0.1 — 0.2 is less than 1 nm. Indeed, for finite 



size ions, w(z) differs from Eq. (36) very close to the in- 
terface and reaches a finite value at z = 0. The size of 
this region exactly corresponds to a which is of the order 
of an ion radius. This is thus an artifact of our point-like 
ion model and occurs only for asymmetric electrolytes at 
neutral surfaces. 

The surface tension a is equal to the excess grand po- 
tential defined as the difference between the grand poten- 
tial of the interfacial system and that of the bulk system: 



parameters (such as a and k^) but also is as close as 
possible to the numerical solution. As suggested by the 
description of Fig.[3j a continuous piecewise <fio(z) is nec- 
essary to account for the essential singularity of p c h(z). 
To show this, let us expand Eq. (41 ) to order O : 



d 2 c 



dz 2 



-4Tr£ B p-q- 
+^-x£BP-q- 



q+e 



e 2 



(43) 



q-e 



This linearization is legitimate, as seen in Fig. [3] 
q+|</>o(z)l < 1 is satisfied for physical valencies. The first 
term in the rhs. of Eq. (43) corresponds to an effective 



local charge source while the second term is responsi- 
ble for the screening of the potential. If we observe the 
charge distribution for q 2 w(z) > q(f>a(z) and z > a, i.e. 



the first term of the rhs. of Eq. (43), we notice that it 



behaves like a distorted peak. The simplest function hav- 



ing a similar behavior is f(z) 



cze 



' K <t> z : where c and 



An 2 


K^ip 2 




32tt 


32ir£ B 










->-] 


>{ 


e 












q+ 







(45) 



-\-iD(z)+q_0o(z) 



The surface tension for electrolytes characterized by q- = 
1 and q + = 1 to 4 is plotted in Fig. [4] function 
of p-, because the anion density is an experimentally 
accessible parameter. Unlike symmetric electrolytes (20) . 
a plot with respect to k 2 may lead to a different behavior. 
One notices that the increase in valency asymmetry leads 
to an important increase of the surface tension. This is of 
course mainly due to the reduction of the cation density 
in the bulk by a factor of q~/q+ necessary to satisfy the 
bulk electroneutrality (see the second term in the integral 
of Eq.pl). 





FIG. 4: (color online) Surface tension d^a/ksT for asymmet- 
ric electrolytes vs. the anion bulk concentration, for increas- 
ing asymmetry q+/q~ = 1 to 4 from bottom to top. 



B. Charged surfaces 

We now consider a symmetric electrolyte in the prox- 
imity of an interface of constant surface charge a s < 
located at z = 0. The variational equation (21 1 simplifies 
to 



d 2 



dz 2 



28{z) 



r^sinh^o- 



(46) 



The mean-field limit (H — > 0) of this equation corre- 
sponds to the NLPB equation, whose solution reads 



<[> (z) = 4arctanh (j b e KbZ ) 



(47) 



where 76 = k& - 



In this Section, we show that a 



piecewise solution for the electrostatic potential similar 



to the one introduced in Section III A agrees very well 
with the numerical solution of Eq. ( 46 1 . Inspired by 



the existence of a salt-free layer close to the interface 
and a mean-field regime far from the interface (WC), 
we propose two types of piecewise variational functions 
(see Appendix [Cj. The first variational choice obeys the 
Poisson equation in the first zone of size h and the non- 
linear Poisson-Boltzmann solution in the second zone : 



4arctanh/7 + 2(z — h) 
4arctanh (j e ~M z -' 1 ) 



for 
for 



z < h, 
z>h, 



(48) 



where 7 = — y 1 + R\ ■ Variational parameters are h 

and an effective inverse screening length k,^. The second 
type of trial potential obeys the Laplace equation with a 
charge renormalization in the first zone and the linearized 
Poisson-Boltzmann solution in the second zone : 



_ 2)7 -k^(z-h) 



for 
for 



z < h, 
5>h. 



(49) 



Variational parameters are h, «0, and the charge renor- 
malization 7/, which takes into account the non-linear ef- 
fects at the mean-field level Q35] ■ The explicit form of the 



FIG. 5: (color online) Electrostatic potential, </>o (in units of 
fcflT): numerical solution of E q. ( |46|l (s ymbols) and restricted 
variational choices Eqs. (48 1 and (I49I) for e = 0, n^c = 4, 



and H = 1, 10, 100, and 1000 (from top to bottom). The vari- 
ational parameters are respectively — 3.83, 3.74, 3.69, 3.66 
and r\ ~ 1. Markers on the rs-axis denote, for each curve, the 
size, h, of the SC zone, plotted vs. InE in the Inset. 



associated variational free energies are reported in Ap- 
pendix [Cj The inset of Fig. [5] displays the size of the SC 
layer h against 3. Our approach predicts a logarithmic 
dependence h oc In S, the factor behind the logarithm 
being k^ 1 for /t& ^> 1. The restricted choices for </>o are 



compared with the full numerical solution of Eq. (41 1 in 
the same figure for e = 0. We see that, as in the pre- 
vious section, the numerical solution and the restricted 
ones match perfectly. Hence salt-exclusion effects are es- 
sentially carried by the parameter h. Furthermore, one 
notices that </>o(-S) relaxes to zero between z = h and 
z = h + 2k~^ > 1 . At k^Iq = 4 we are in the linear regime 
of the PB equation and therefore one has r\ ~ 1. The 
charge renormalization idea was introduced by Alexan- 
der et al. [14], who showed that the non-linearity of the 
PB equation can be effectively taken into account at long 
distances by renormalizing the fixed charge source and 
extending the linearized zone where \(f>o\ < 1 to the whole 
domain. A linear solution of the form Eq. ( [49] ) can be 
very helpful for complicated geometries or in the presence 
of a non-uniform charge distribution where the NLPB 
equation does not present an analytical solution even at 
the mean-field level. These issues will be discussed in a 
future work. 

Figure [6 displays the ion concentrations Pi (z)l Pi.b = 
e - **, which are related to the ion PMF Eq. (24), com- 



puted with the restricted solution Eq. ( 48 1 for several val 



ues of S. As already said in the Introduction, in rescaled 
distance, the coupling parameter S measures the strength 
of the excess chemical potential, w(z). We first see that 
for coions as well as for counterions, the depletion layer 
in rescaled units in the proximity of the dielectric inter- 
face increases with S due to the image charge repulsion 
and/or solvation effect, i.e. the term e~^ w ^ in Eq. (46 1. 



Furthermore, one notices that the counterion density ex- 
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(a) 



P'Pb 




(b) 

FIG. 6: (color online) Ion densities for KblG = 4, and (a) 
e = and (b) e = e w , for increasing coupling parameter: from 
left to right, S = 1, 10, 100, and 1000. Solid lines correspond 
to counterions, dashed lines to coions and dashed-dotted lines 
to the Poisson-Boltzmann result (47 1. 



hibits a maximum. This concentration peak is due to 
the competition between the attractive force towards the 
charged wall and the repulsive image and solvation in- 
teractions. It is important to note that in the particular 
case e = e w , there is no depletion layer for S < 10. 



IV. DOUBLE INTERFACE 

In this Section, the variational method is applied to a 
double interface system which consists of a slit-like pore 
of thickness d, in contact with an external ion reservoir 
at its extremities (Fig. [lj . The dielectric constant is e w 
inside the pore and e in the outer space. The electrolyte 
occupies the pore and the external space is salt-free. The 



solution of the DH equation ( A2 ) in this geometry is [6] 



w(z) = (n b - K v )£l 



(50) 



+ Ib 



kdk 



A(k/K v ) 



yFT^ e 2V^^_A2(fc/ Kll ) 



x 2A(fc/ K ,) + e 2 ^ 2 )^^? +e 2z^+^ 



where A(x) is given in Eq. (33). The variational pa- 
rameter of the Green's function is the variational inverse 
screening length k v which is taken uniform (generalized 



Onsager-Samaras approximation, see [S1[5T]). A more 
complicated approach has been previously developed in 
Ref. [3T] where the authors introduced a piecewise form 
for the variational screening length, i.e. k(z) — k v over 
a layer of size h and n v — Kb in the middle of the pore. 
Although this choice is more general than ours, the min- 
imization procedure with respect to k v is significantly 
longer than in our case and the variational equation is 
much more complicated. Consequently, this piecewise 
approach is not very practical when one wishes to study 
a charged membrane where the external field created by 
the surface charge considerably complicates the technical 
task (see Section IV B[ ). We show that the simple varia- 
tional choice adopted here captures the essential physics 
with less computational effort. 

As in Eq. (32), the integral on the rhs. of Eq. (50) 



takes into account both image charge and solvation ef- 
fects due to the two interfaces, whereas the first term 
is the Debye result for the difference between the bulk 
and a hypothetic bulk of inverse screening length k v . We 
should emphasize that, in the present case, the spatial in- 
tegrations in Eqs. (A3)-(A4) run over the confined space, 
that is from z — to z = d. By substituting the so- 
lution Eq. (50) into Eqs. p0|-( A5 1 and performing the 
integration over z, one finds [22] 



24tt 



16tt 

dxxln [l - A 2 (x)e~ 2 ^ dx ] 



F 2 + F 3 
S 

K v 

47T . n 

1 f°° (A(x)-A 3 (x)) /x -2n v dA 2 (x) 

L dx 



(51) 



~8ir 



e 



2ih, 



where we have defined A(x) = A (\fx 2 — l). 

The limiting case e = allows for closed-form expres- 
sions. This limit is a good approximation for describing 
biological and artificial pores characterized by an exter- 
nal dielectric constant much lower than the internal one. 
In the following part of the work, we will deal most of the 
time with the special case e = 0, unless stated otherwise. 
In this limit, Eq. (51) simplifies to 



Fo + Fa 



KyCl 

24^ 
8^ 



16tt L 



2 In ( 1 



Li 2 (e- 2d ^) 



Li 3 (e 



-2dK„ 



)] 



167rd 2 



(52) 



where Li n (x) stands for the polylogarithm function and 
£(x) the Riemann zeta function (see Appendix |d| . 
Within the same limit (e = 0), A(x) = 1 and we obtain 
an analytical expression for the Green's function Eq. ( 50 ) 



w (z) = (Kb — k v )£b -T- In (l 

+ ^ 
2d 



-2dn 



pf e -2d Kv . 



d Our , „ / Z 

-e- 2 ^ z 2 F 1 1,-,1 
z V d 



-2dK 



(53) 
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where f3(x;y,z) is the incomplete Beta function and 
2Fi(a, b; c; d) the hypergeometric series. The definitions 
of these special functions are given in Appendix [D] At 
this step, the PMF thus depends on three adimensional 
parameters, namely dn v , dnt, and d/is- 

For the system with a single interface, the ion fugacity 
Xi was fixed by the bulk density. In the present case 
where the confined system is in contact with an external 
reservoir, Ai is fixed by chemical equilibrium: 



A* = Xi,b = Pi.be- 2 



(54) 



where Kb and Xt t b are respectively the inverse Debye 
screening length and the fugacity in the bulk reservoir 
[see Eq. (28)]. Once this constraint is taken into ac- 
count, the last term of electrostatic part of the vari- 
ational grand potential Eq. (20 1 can be written as 

Eq. (21) then becomes for a symmetric q : q elec- 
trolyte: 



d 2 



dz 2 

The optimization of F v 



W*) sinh0 o = -4nrqi B tT. [S(z) + 5{z - d)] 

(551 

Ft + F 2 + F 3 given by Eq. ( 20 ) 



and (52 ) with respect to the inverse trial screening length 



k v leads to the following variational equation for k v : 



(d Kv ) 2 + dK v tar>h(dK v ) = (dK b ) 2 / dxe-^ w ^ xd) 



x co S h^M](l + ^|^ | -(56) 
|_ cosh(dK„) J 

Within the particular choice that fixed the functional 



form of the k v dependent Green's function Eq. (53), the 



two coupled equations (55 1 and (|56[) are the most general 



variational equations. In the following, we first consider 
the case of neutral pores and then the more general case 
of charged pores. 



A. Neutral pore, symmetric electrolyte 

In the case of a symmetric q : q electrolyte and a neu- 
tral membrane, a s — 0, the solution of Eq. ( 55 ) is natu- 
rally 0o = 0. The variational parameter solution of 
Eq. (56) with fa — and w(z) — wq(z) given by Eq. (53) 
when e = 0, which can be written as dn v = f(dKh, ^b/HJ. 



Let us note that Eq. ( 56 1 can be solved with the Mathe- 



matica software in a fraction of a second. 

Within the Debye-Hiickel closure approach, 
Yaroshchuk (see Eq. (59) of Ref. [5]) obtains a self- 
consistent approximation for constant k v by replacing 
the exponential term of Eq. (12) with its average value 
in the pore: 



da; ( 



'j(xd) 



(57) 




FIG. 7: (color online) Inverse screening length inside the neu- 
tral membrane (monovalent ions) normalized by Kb vs. the 
pore size d/ls for e = and (a) k^Ib = 0.1 (pb = 1.926 mM), 
(b) KblB = 1 {pb = 0.1926 M). Dashed lines correspond to 
the mid-point approximation, Eq. (58 1. The inset shows the 



characteristic pore size corresponding to total ionic exclusion 
as a function of the inverse bulk screening length. The bot- 
tom curve corresponds to monovalent ions and the top curve 
to divalent ions. 



which should be compared with Eq. ( 56 ) with fa = 0. In 
order to simplify the numerical task, Yaroshchuk intro- 
duces a further approximation in which he replaces the 
potential w(z) inside the depletion term of Eq. (57) by its 
value in the middle of the pore, w(d/2). Then Eq. (57) 
takes the simpler form 



u(d/2) 



(58) 



The self-consistent midpoint approximation is frequently 
used in nanofiltration theories [51 [55J [35] . For e = 0, the 
mid-point potential has the simple form w(d/2) = (kj — 
k v )£b — 2£b hi(l — e~ Kvd )/d. This approach is compared 
with the full variational treatment in Fig. [7] where the 
adimensional inverse screening length in the pore k v / kj, 
is plotted as a function of the pore size d. We first note 
that as d decreases below a critical value d* , the pore is 
empty of salt and k v = 0. The inset of Fig. [7] shows d* 
versus the inverse bulk screening length. Searching for d 
such that k v = in Eq. ( |56[ ) leads to the same equation 
as Eq. (571, thus the value of d* is identical within both 
approaches. However, Fig. [7] shows that the mid-point 
approximation, Eq. ( 58 ) , overestimates the internal salt 



concentration as well as the abruptness of the crossover 
to an ion-free regime for decreasing pore size. Indeed, 
this approximation is equivalent to neglecting the strong 
ion exclusion close to the pore surfaces (which is larger 
than in the middle of the pore). A similar behavior was 
also observed in Fig. 6 of Ref. [5T] for the screening length 
in the neighborhood of the dielectric interface. 

The effect of the dielectric discontinuity is illustrated 
in Fig. ^fa) where the inverse internal screening length is 
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(a) 





(b) 
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FIG. 8: (color online) Inverse screening length inside the 
membrane vs. the pore size cL/Ib (ew = 78, k^Ib = 1). (a) 
From bottom to top: e = (A = 1), e = 3.2 (A = 0.92), 
e = 39 (A = 1/3), and e = 78 (A = 0). (b) Log-linear plot 
for monovalent, divalent and trivalent ions, from left to right 
(e = 0). 



compared for e between and e w — 78 where the image- 
charge repulsion is absent and the solvation effect is solely 
responsible for ion repulsion. First of all, one observes 
that the total exclusion of ions in small pores is specific to 
the case e = 0. Moreover, in the solvation only case, the 
inverse screening length inside the pore only slightly de- 
viates from the bulk value, 0.8 < K v /Kb < 1. This clearly 
indicates that, within the point-like ion model consid- 
ered in this work, the image-charge interaction brings the 
main contribution to salt-rejection from neutral mem- 
branes. Roughly speaking, the image-charge and solva- 
tion effects come into play when the surface of the ionic 
cloud of radius n^ 1 around a single ion located at the pore 
center touches the pore wall, i.e. for k^ 1 > d/2. This 
simple picture fixes a characteristic length d c h — 2k^~ 1 
below which the internal ion density significantly devi- 
ates from the bulk value and ion-rejection takes place. 
This can be verified for intermediate salt densities in the 
bottom plot of Fig. [7] and the top plot of Fig. [8] 

Since image-charge effects are proportional to q 2 , we 
illustrate in Fig. Mb) the effect of ion valency q. At 
pore size d ~ 2.5£g ~ 1.8 nm, where the inverse inter- 
nal screening length for monovalent ions is close to 80 % 
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FIG. 9: (color online) Salt reflection coefficient (dimension- 
less) against the logarithm of the inverse bulk screening length 
for e = and two pore sizes, djis =2 and 5 (red lines corre- 
spond to the mid-point approximation, Eq. (58 1). 



of its saturation value kj,, the exclusion of divalent ions 
from the membrane is total. This effect driven by image 
interactions is even much more pronounced for trivalent 
ions. Since the typical pore size of nano-filtration mem- 
branes ranges between 0.5 and 2 nm, we thus explain why 
ion valency can play a central role in ion selectivity, even 
inside neutral pores. 

The salt reflection coefficient, frequently used in mem- 
brane transport theories to characterize the maximum 
salt-rejection (obtained at high pressure) is related to the 
ratio of the net flux of ions across the membrane to that 
of the solvent volume flux J per unit transverse surface : 



1 



1 

JPb Jo 



v\\(z)p(z)dz 



(59) 



= 1-12 / x(l-x) e- q - w{xd) dx. 

Jl/2 

where we have used, in the second equality, the Poiseuille 
velocity profile, vn(z) — ^p-z(d — z) in the pore and the 
PMF given by Eq. (24 1. It depends only on the parame- 
ters Kb(-B and d/£s- In certain nanopores with hydropho- 
bic surfaces, the solvent flux may considerably deviate 
from the Poiseuille profile (see [37]). In this case, the ve- 
locity profile is flat, Vn(z) = 4. We emphasize that since 
the velocity profile is normalized in both cases, the mid- 
point approximation is unable to distinguish between a 
Poiseuille and a plug flow velocity profile. Fig. [9] displays 
S s as a function of the inverse bulk screening length 
for two pore sizes d = 2£b and d = 5£b- As seen by 
Yaroshchuk, decreasing the pore size shifts the curves to 
higher bulk concentration and thus increases the range of 
bulk concentration where nearly total salt rejection oc- 
curs. However, quantitatively, the difference between the 
variational and mid-point approaches becomes significant 
at high bulk concentrations and this difference is accen- 
tuated in the case of plug-flow (for which E s is higher 
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when compared to the Poiseuille case because the flow 
velocity no longer vanishes at the pore wall where the 
salt exclusion is strongest). This deviation is again due 
to the midpoint approximation of Eq. ( 58 1 in which the 
image interactions are underestimated. However since 
the velocity profile vanishes at the solid surface for the 
Poiseuille flow, the deficiencies of the mid-point approx- 
imation are less visible in E s than in k v in this case. 

Finally, we compute the disjoining pressure within our 
variational approach. We compare in Appendix [E] the 
result with that of the more involved variational scheme 
presented in Ref. [H] and show that one gets a very 
similar behaviour, revealing that the simpler variational 
method is able to capture the essential physics of the slit 
pore. 

As stressed above, the main benefit obtained from the 
simpler approach proposed in this work is that the min- 
imization procedure is much less time consuming. This 
point becomes crucial when considering the fixed charge 
of the membrane, which is thoroughly studied in the next 
section. 
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FIG. 10: (color online) Inverse internal screening length k v 
against 3 for Kb£a — 2, e = and (a) d = 3£g and (b) d = 
10£g- Comparison of various approximations: Yaroshchuk, 



Eq. (71 1 (diamonds), variational Donnan potential (dashed 
line), piecewise solutions (solid line), and numerical results 
(squares). Horizontal lines corresponds to the WC limit, 



Eq. 1671 (top), and SC limit, Eq. (70 1 (bottom) 



B. Charged pore, symmetric electrolyte 

In this section, we apply the variational approach to a 
slit-like pore of surface charge a s < 0. In the following, 
we will solve Eqs. (55) and (56) numerically in order to 



test, as in the case of a single charged surface, the validity 
of restricted trial forms for (f>o(z). We define the partition 
coefficients in the pore for countcrions and coions, k+ and 
k , as 



k+ = 



P± 
Pb 



dz 



-*±(*) 



(60) 



where <&±(z) is given by Eq. (26) 



1. Effective Donnan Potential 

When one considers a charged nanopore, because of its 
small size, gradients of the potential (j>o can be neglected 
as a first approximation. We thus assume a constant po- 
tential (pQ. The so-called effective Donnan potential 4> 
introduced by Yaroshchuk [5] will be fixed by the vari- 
atio nal p rinciple. By differentiating the grand potential 
Eq. (|20| ) with respect to <f>o (or equivalently integrating 
Eq. (|55[) from z = to z = d with V0o = 0), we find 



2 |<r s | = -2qp b smh(qct> ) I dz e 2 
Jo 



(61) 



which is simply the electroneutrality relation in the pore, 
taken in charge by the electrostatic potential 4>q. By 



defining 



/ dxexp[- 
Jo 



■q 2 £ B w(xd)/(2d)} 



dx exp[— , Ew(xd) / (2d)], 



(62) 
(63) 



where w(x) = w{x)d/tB, we have k± — T exp(=Fg^o) and 
Eq. (61) can be rewritten as 



qpbd 



x m 
qpb 



n 2 h dl G 



k 2 b d 



(64) 



where the second equality contains the Gouy-Chapman 
length Iq and the quantity X m — 2\o~ s \/d, frequently 
used in nanofiltration theories, corresponds to the volume 
charge density of the membrane. Hence, the partition 
coefficient of the charge, k + — does not depend on E, 



i.e. charge image and solvation forces. By using Eq. (61 ) 
in order to eliminate the potential (f>o from Eq. (60), one 
can rewrite the partition coefficients in the form 



k+ = Te Tq t° = 



\ 



r 2 



± 



(65) 



By substituting into Eq. fl56| ) the ana lyti cal expression 
for <f>o obtained from Eq. (61) (or Eq. (65)), one obtains 
a single variational equation for k v to be solved numeri- 
cally, 



(dk v ) 2 + dk v ta,nh(dk v ) 



dx e 




rcosh(c?K„) 
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FIG. 11: (color online) Ionic partition coefficients, k±, vs. S 
for K b £ G = 2, e = 0, and (a) d = 3£ G and (b) d = IQ£ G - The 
horizontal line corresponds to the SC limit for counterions. 
As explained in the text, we note that fc+ — fc_ = 8/(k%cI£q). 



The numerical solution of Eq. (66 1 is plotted in Fig. 10 as 



a function of the coupling parameter S. We see that as we 
move from the WC limit to the SC one by increasing S, 
the pore evolves from a high to a low salt regime. This 
quite rapid crossover, which results from the exclusion 
of ions from the membrane, is mainly due to repulsive 
image-charge and solvation forces controlled by T whose 
effects increase with increasing 5. 

In Fig. [TT] are plotted the partition coefficients of coun- 
terions and coions, Eq. ( [60] ), as a function of S. Here 
again, k± decreases with increasing S. Moreover, we 
clearly see that the rejection of coions from the mem- 
brane becomes total for 2 > 4. In other words, even 
for intermediate coupling parameter values, we are in a 
counterion-only state. This is obviously related to the 
electrical repulsion of coions by the charged surface. 

In the asymptotic WC limit (S 0), T = 1 and we find 
the classical Donnan results in mean-field where fe_ = 
k^ 1 = e q ^" with q<f> Q — arcsinh[4/(/«j*d)]. The variational 



equations ( 66 ) and ( 65 ) reduce to 



k+ = 



\ 



± 



(67) 
(68) 



Quite interestingly, the relation Eq. (67 1 shows that, even 



in the mean-field limit, due to the ion charge imbalance 
created by the pore surface charge, the inverse screening 
length is larger than the Debye-Hiickel value k&. In the 
case of small pores or strongly charged pores or at low 
values of the bulk ionic strength, i.e. K%lQd < 1 or 
dpi, <C \cr s \/q, we find k v ~ 2/y/Icd and p_ = and p+ = 



2\a s \/(dq). We thus find the classical Poisson-Boltzmann 
result for counterions only |24) . The counterion-only case 
is also called good coion exclusion limit (GCE), a notion 
introduced in the context of nanofiltration theories [51 
|3"51 [313] . Hence, in this limit the quantity of counterions 
in the membrane is independent of the bulk density and 
depends only on the pore size d and the surface charge 
density o~ s . In the case of a membrane of size d ~ 1 nm 
and fixed surface charge a s ~ 0.03 nm~ 2 , this limit can 
be reached with an electrolyte of bulk concentration pb ~ 
50 mM. In the opposite limit K^lad ^> 1, one finds 
k v ~ Kb and p± — pb- 

In the SC limit S — >■ oo, r = and Eq. (66 1 simplifies 

to 

(dk v ) 2 + dk v tanh(dK„) = Ad[l + sech.(dh v )]. (69) 



For d > £g (d > 1), the solution of Eq. (69 1 yields with 
a high accuracy 



Vl + 16 J - 1 



2d 



(70) 



The partition coefficients simplify to fc_ = and k + — 
8/(dkl) = 2\a s \/(dqpi,) and we find the counterion only 
case (or GCE limit) without image charge forces dis- 
cussed by Netz [21] . Partition coefficients in the SC limit 
and variational inverse screening length in both limits, 
Eqs. (67) and (70 1, are illustrated in Figs. 10 and 
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by dotted reference lines. Consequently, one reaches for 
H = the GCE limit exclusively for low salt density 
or small pore size, while the SC limit leads to GCE for 
arbitrary bulk density. It is also important to note that 
although the pore-averaged densities of ions are the same 
in the GCE limit of WC and SC regimes, the density pro- 
files are different since when one moves away from the 
pore center, the counterion densities close to the inter- 
face increase in the WC limit due to the surface charge 
attraction and decrease in the SC limit due to the image 
charge repulsion. 

It is interesting to compare this variational approach to 
the approximate mid-point approach of Yaroshchuk [B]. 
For charged membranes, he considers a constant po- 
tential and replaces the exponential term of Eqs. (11) 



and ( 12 ) by its value in the middle of the pore. He ob- 



tains the following self-consistent equations: 



K 

2 k,| 



-2qdp b sinh (qfa) g-^W 2 ). 



(71) 
(72) 



The above set of equations are frequently used in nanofil- 
tration theories [QEISJES]- By combining these equations 
in order to eliminate 4>o, one obtains an approximate non- 
linear equation for k v (approximation CYar in Fig. 10). 



In the limit of a high surface charge, the non-linear equa- 
tions (l7Tj)-([72|) depend only on the pore size d and the 



surface charge density a s 



8ir£ B q\o- s 



4 



(73) 
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FIG. 12: (color online) Inverse internal screening length k v 
against the reduced surface charge a = £ B a s for d = £b, e = 
and (a) k\>£b = 1, (b) Kbla = 2: constant variational Donnan 



approximation (solid line), asymptotic result Eq. (731 (dotted 
line) and Yaroshchuk approximation Eq. (71 1 (dashed line). 



One can verify that in the regime of strong surface 
charge, Eq. (73) is also obtained from the asymptotic 
solution Eq. ( 70 ) since the dependence of the PMF on 
z is killed when 3 —¥ oo and only the mid-pore value 
contributes. The numerical solution of Eq. ( 71 1-( 72 ) 
is illustrated as a function of E in Fig. [TUJ and as a 
function of the surface charge in Fig. [l2j together with 
the asymptotic formula Eq. ( 73 1 . For the parameter- 
range considered in Fig 



10 the solution of Eq. (71) 



strongly deviates from the result of the full variational 
calculation. For S < 2, the mid-point approach follows 
an incorrect trend with increasing S. It is clearly 
seen that at some values of the coupling parameter, 
Eqs. (71)-(72) do not even present a numerical solution. 
Using the relations cL/£b = d/S and f^Kj = for 
monovalent ions, one can verify that the regime where 
the important deviations take place corresponds to high 
ion concentrations. This is confirmed in Fig. 12 the 
error incurred by the approximate mid-point solution of 
Yaroshchuk increases with the electrolyte concentration. 



In Section [IV A| on neutral nanopores, it has been un- 
derlined that, due to the image charge repulsion, the ionic 
concentration inside the pore increases with the pore size 
d (see Fig. [8]). In the present case of charged nanopores, 
this result is modified: Eqs. ( 67 ) , ( 70 ) and ( 73 1 show that 



for strongly charged nanopores the concentration of ions 
inside the pore decreases with d. Moreover, the very high 
charge limit is a counterion-only state and Eq. (61 ) shows 



that, for a fixed surface charge density, electroneutrality 
alone fixes the number of counterions, AL., in a layer 
of length d joining both interfaces, and image charge 
interactions play a little role. This is the reason why 
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FIG. 13: (color online) Partition coefficient in the pore 
of coions (a) and counterions (b) vs. the pore size 
d/£ B for increasing surface charge density, a B £ B = 
0, 0.004, 0.02, 0.04, 0.08, 0.12, from left to right, and n b £ B = 1. 
Inset: Critical pore size d CI vs. the surface charge density a s 
(e = 0). 



Hence, we expect an intermediate charge regime which 
interpolate between image force counterion repulsion 
(case of neutral pores, see Section IV A) and counterion 
attraction by the fixed surface charge. This is illustrated 
in Fig. [T3| where the partition coefficients are plotted vs. 
d for increasing a s . As expected, coions are electrostat- 
ically pushed away by the surface charge which adds to 
the repulsive image forces, leading to a stronger coion ex- 
clusion than for neutral pores. The issue is more subtle 
for counterions: obviously, increasing the surface charge, 
o s , at constant pore size, d, increases However, for 
small fixed a Sl a regime where image charge and direct 
electrostatic forces compete, k + is non-monotonic with 
d. Below a characteristic pore size, d < d CI , the electro- 
static attraction dominates over image charge repulsion 



K v OC P- 



N-/(Sd) decreases for increasing d. 



and due to the mechanism explained above, fc+ decreases 
for increasing d. For d > d CI , the effect of the surface 
charge weakens and k + starts increasing with d. In this 
regime, the pore behaves like a neutral system. The inset 
of Fig.[T3]shows that d CI increases when a s increases. For 
highly charged membranes l 2 B a s ^> 0.1, there is no min- 
imum in k+(d), and the average counterion density in- 
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FIG. 14: (color online) Ion densities in the nanopore for e = 
ew, 3 = 1 and H/Ig = 2. The continuous lines correspond to 
the prediction of the variational method with four parameters, 
the dashed-dotted line the variational solution with a single 
parameter (see the text), the symbols are MC results (Fig. 2 
of [29]) and the dashed lines denote the numerical solution of 
the non-linear PB result. 



side the membrane monotonically decreases towards the 
bulk value. Experimental values for surface charges are 
< <t s < 0.5 nm~ 2 (or < £%a s < 0.25), which cor- 
responds to physically attainable values of d cr . The in- 
terplay between image forces and direct electrostatic at- 
traction is thus relevant to the experimental situation. 

The variational Donnan potential approximation is 
thus of great interest since it yields physical insight into 
the exclusion mechanism and allows a reduction of the 
computational complexity. However, membranes and 
nanopores are often highly charged and spatial variations 
of the electrostatic potential inside the pore may play an 
important role. In the following we seek a piecewise so- 
lution for 0o ( z )- 



2. Piecewise solution 



The variational modified PB equation ( 55 1 for 0o shows 
that as one goes closer to the dielectric interface, w(z) 
increases and the screening experienced by the poten- 
tial </>o gradually decreases because of ionic exclusion. 
This non-perturbative effect which originates from the 
strong charge-image repulsion inspires our choice for the 
variational potential 4>o{z). We opt for a piecewise so- 
lution as in Section Hilt a salt-free solution in the zone 
< z < h and the solution of the linearized PB equa- 
tion for h < z < d/2, with a charge renormalization 
parameter 77 taking into account non-linear effects. By 
inserting the boundary conditions dtpa/dz\ z= o = 2t)I£q 
and d4>o/dz\ z=d / 2 = and imposing the continuity of </>o 
and its first derivative at z = h [Eq. (B3)], the piecewise 
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FIG. 15: (color online) Variational electrostatic potential (in 
units of feT) in the nanopore. Comparison of the numeri- 
cal solution of Eq. ( |55[ l with the homogeneous (h = 0) and 
piecewise solution of Eq. (C5l for icn-b = 3 and (a) 3 = 1, 
(b) 3 = 100. The ho rizontal line is the Donnan potential 
obtained from Eq. (61 1 (e = 0). 



potential, solution of Eq. (55) with ft? exp[— q 2 w(z)/2] 



replaced by k£, takes the form 



<f>o( z ) 
where 



for < z < h, 



in- ?a U_ d\ 
V la\ 2 I 

w ?H- cosh[K <t ,(d/2-z)] r , < < , 



coth 



/. , . ( ^ h 



(74) 



(75) 



is imposed by continuity, and kq, <p, h and r\ are the vari- 
ational parameters. By injecting the piecewise solution 
Eq. (|C5|) into Eq. pOl, we finally obtain 



S 



2K 

q 

V(V " 4 ) 

2£ G K c j ! 

A-K£ B q 2 



77(77-2) 



h 



V 2 (d/2 - h) 



l G 2£ G smh 2 [n 4> (d/2- h)} 
coth[K (d/2 - h)] + (ysj 

(76) 



dz e'^ w( - z) cosh0 o (z). 



The solution to the variational problem is found by op- 
timization of the total grand potential F = F\ + F2 with 
respect to k^, ip, h, 77 and k v , where F2 + F3 is given by 
Eq. ( 52 ) for a general value of e and by Eq. ( 52 ) for e = 0. 



This was easily carried out with Mathematica software. 

A posteriori, we checked that two restricted forms for 
(f>o, homogeneous with h = and piecewise with ip = 0, 
were good variational choices. Fig. [M] compares the ion 
densities obtained from the variational approach (with 
homogeneous </>o) with the predictions of the MC simu- 
lations [33] and the NLPB equation for e = e^, d = 2 
and 3 = 1. Two variational choices are displayed in this 
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figure, namely, the homogeneous approach with four pa- 
rameters k v — 1.68, k^, = 1.36, ip — 0.16, 77 = 0.97 and a 
simpler choice with rj = 1, k$ = K v and two variational 
parameters: k v — 1.69,0 = —0.18. In the latter case, 
one can obtain an analytical solution for tp and injecting 
this solution into the free energy, one is left with a single 
parameter n v to be varied in order to find the optimal 
solution. We notice that with both choices, the agree- 
ment between the variational method and MC result is 
good. It is clearly seen that the proposed approach can 
reproduce with a good quantitative accuracy the reduced 
solvation induced ionic exclusion, an effect absent at the 
mean-field level. Moreover, we verified that with the sin- 
gle parameter choice, one can reproduce at the mean-field 
variational level the ion density profiles obtained from 
the numerical solution of the NLPB equation (dashed 
lines in Fig. 14 ) almost exactly. We finally note that the 



small discrepancy between the predictions of the varia- 
tional approach and the MC results close to the interface 
may be due to cither numerical errors in the simulation, 
or our use of the generalized Onsager-Samaras approx- 
imation (our homogeneous choice for the inverse effec- 
tive screening length appearing in the Green's function 
vq does not account for local enhancement or diminution 
of ionic screening due to variations in local ionic density) . 

For e = 0, the piecewise and homogeneous solu- 
tions are compared with the full numerical solution of 
Eqs. (|55}-(|56) in Fig. pj for 3 = 1 and 100. First of 
all, one observes that for 3 = 1, both variational solu- 
tions match perfectly well with the numerical solutions. 
For 3 = 100, the piecewise solution matches also per- 
fectly well with the numerical one, whereas the match- 
ing of the homogeneous one is poorer. The optimal val- 
ues of the variational parameters (k v , K<f,,r],h) for the 
piecewise choice are (2.57,2.6,0.98,0.15) for 3=1 and 
(0.83, 0.13, 0.97, 1.37) for S = 100. 

The form of the electrostatic potential <fio(z) is inti- 
mately related to ionic concentrations. Ion densities in- 



side the pore are plotted in Fig. 16 for 3 = 1 and 2 = 100. 
We first notice that even at 3 = 1, the counterion density 
is quite different from the mean-field prediction. Further- 
more, due to image charge and electrostatic repulsions 
from both sides, the coion density has its maximum in 
the middle of the pore. On the other hand, the coun- 
terion density exhibits a double peak, symmetric with 
respect to the middle of the pore, which originates from 
the attractive force created by the fixed charge and the 
repulsive image forces. When 3 increases, we see that the 
counterion density close to the wall shrinks and becomes 
practically flat in the middle of the pore. Hence the po- 
tential (f>o linearly increases with z until the counterion- 
peak is reached and then it remains almost constant since 
the counterion layer screens the electrostatic field created 
by the surface charge (since in Fig. 15 z is renormalized 



by the Gouy-Chapman length which decreases with in- 
creasing o s , one does not see the increase of the slope 
4>o(z — 0)). In agreement with the variational Donnan 
approximation above, coions are totally excluded from 
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FIG. 16: (color online) Local ionic partition coefficient in the 
nanopore (same parameters as in Fig. 15 and e = 0) com- 
puted with the piecewise solution, (a) H = 1 , (b) S = 100. 
The dotted line in the top plot corresponds to the mean-field 
prediction for counterion density. 



the pore for large 3. Hence, the piecewise potential al- 
lows one to go beyond the variational Donnan approxi- 
mation within which the density profile does not exhibit 
any concentration peak. 

The inverse screening length k v obtained with the 



piecewise solution is compared in Fig. 10 with the pre- 
diction of the Donnan approximation and that of the 
numerical solution. The agreement between piecewise 
and numerical solutions is extremely good. Although 
the Donnan approximation slightly underestimates the 
salt density in the pore, its predictions follow the correct 
trend. 



V. CONCLUSION 

In this study, we applied the variational method to 
interacting point-like ions in the presence of dielectric 
discontinuities and charged boundaries. This approach 
interpolates between the WC limit (S < 1) and the SC 
one (3 3> 1), originally defined for charged boundaries 
without dielectric discontinuity, and takes into account 
image charge repulsion and solvation effects. The vari- 
ational Greens 's function vo has a Debye-Hiickel form 
with a variational parameter k v and the average vari- 
ational electrostatic potential 4>q{z) is either computed 
numerically or a restricted form is chosen with varia- 
tional parameters. The physical content of our restricted 
variational choices can be ascertained by inspecting the 
general variational equations (Eqs. ( 11 1-( 12 1 for symmet- 
ric salts). The generalized Onsager-Samaras approxi- 
mation that we have adopted for the Green's function 
replaces a local spatially varying screening length by a 
constant variational one; although near a single interface 
this screening length is equal to the bulk one, in confined 
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geometries the constant variational screening length can 
account in an average way for the modified ionic environ- 
ment (as compared with the external bulk with which the 
pore is in equilibrium) and can therefore strongly devi- 
ate from the bulk value. This modified ionic environ- 
ment arises both from dielectric and reduced solvation 
effects present even near neutral surfaces (encoded in the 
Green's function) and the surface charge effects encoded 
in the average electrostatic potential. Our restricted vari- 
ational choice for (j>o(z) is based on the usual non- linear 
Poisson-Boltzmann type solutions with a renormalized 
inverse screening length that may differ from the one 
used for vq and a renormalized external charge source. 
The coupling between Vq and </>o arises because the in- 
verse screening length for vq depends on <fi and vice- 
versa. The optimal choices are the ones that extremize 
the variational free energy. 

In the first part of the work, we considered single in- 
terface systems. For asymmetric electrolytes at a single 
neutral interface, the potential <j>o(z) created by charge 
separation was numerically computed. It was satisfac- 
torily compared to a restricted piecewise variational so- 
lution and both charge densities and surface tension are 
calculated in a simpler way than Bravina [5] and valid 
over a larger bulk concentration range. The variational 
approach was then applied to a single charged surface and 
it was shown that a piecewise solution, characterized by 
two zones, can accurately reproduce the correlations and 
non-linear effects embodied in the more general varia- 
tional equation. The first zone of size h is governed by a 
salt-free regime, while the second region corresponds to 
an effective mean-field limit. The variational calculation 
predicts a relation between h and the surface charge of 
the form h oc (c + In |er s |)/|<7 s | where the parameter c 
depends on the temperature and ion valency. 

In the second part, we dealt with a symmetric elec- 
trolyte confined between two dielectric interfaces and in- 
vestigated the important problem of ion rejection from 
neutral and charged membranes. We illustrated the ef- 
fects of ion valency and dielectric discontinuity on the 
ion rejection mechanism by focusing on ion partition and 
salt reflection coefficients. We computed within a vari- 
ational Donnan potential approximation, the inverse in- 
ternal screening length and ion partition coefficients, and 
showed that for S > 4 one reaches the SC limit, where 
the partition coefficients are independent of the bulk con- 
centration and depend only on the size and charge of the 
nanopore. This result has important experimental appli- 
cations, since it indicates that complete filtration can be 
done at low bulk salt concentration and/or high surface 
charge. Furthermore, we showed that, due to image inter- 
actions, the quantity of salt allowed to penetrate inside a 
neutral nanopore increases with the pore-size. In the case 
of strongly charged membranes, this behavior is reversed 
for the whole physical range of pore size. We quanti- 
fied the interplay between the image charge repulsion and 
the surface charge attraction for counterions and found 
that even in the presence of a weak surface charge, the 



competition between them leads to a characteristic pore 
size d cr below which the counterion partition coefficient 
rapidly decreases with increasing pore size. On the other 
hand, for nanopores of size larger than d cr the system be- 
haves like a neutral pore. Our variational calculation was 
compared to the Debye closure approach and the mid- 
point approximation used by Yaroshchuk 6J. The clo- 
sure equations have no exact solution even at the numer- 
ical level. Our approach, based on restricted variational 
choices, shows significant deviations from Yaroshchuk's 
mid-point approach at high ion concentrations and small 
pore size. Finally, the introduction of a simple piecewise 
trial potential for <f>Q, which perfectly matches the nu- 
merical solutions of the variational equations, enabled us 
to go beyond the variational Donnan potential approxi- 
mation and thus account for the concentration peaks in 
counterion densities. We computed ion densities in the 
pore and showed that for S > 4, the exclusion of coions 
from the pore is total. We also compared the ionic den- 
sity profiles obtained from the variational method with 
MC simulation results and showed that the agreement is 
quite good, which illustrates the accuracy of the varia- 
tional approach in handling the correlation effects absent 
at the mean-field level. 

The main goal in this work was first to connect two 
different fields in the chemical physics of ionic solutions 
focusing on complex interactions with surfaces: field- 
theoretic calculations and nanofiltration studies. More- 
over, on the one hand, this variational method allows one 
to consider, in a non-perturbative way, correlations and 
non-linear effects; on the other hand the choice of one 
constant variational Debye-Hiickel parameter is simple 
enough to reproduce previous results and to illuminate 
the mechanisms at play. This approach is also able to 
handle, in a very near future, more complicated geome- 
tries, such as cylindrical nanopores, or a non-uniform sur- 
face charge distribution. 

The present variational scheme also neglects ion-size 
effects and gives rise to an instability of the free energy 
at extremely high salt concentration. Second order 
corrections to the variational method may be necessary 
in order to properly consider ionic correlations leading 
to pairing [3U1 HI] and to describe the physics of charged 
liquids at high valency, high concentrations or low 
temperatures. Introducing ion size will also allow us to 
introduce an effective dielectric permittivity e p for water 
confined in a nanopore intermediate between that of 
the membrane matrix and bulk water, leading naturally 
to a Born-self energy term that varies inversely with 
ion size and depends on the difference between l/e w 
and l/e p [421 143] . Furthermore, the incorporation of 
the ion polarizability |44j will yield a more complete 
physical description of the behavior of large ions [45] . 
Charge inversion phenomena for planar and curved 
interfaces is another important phenomenon that we 
would like to consider in the future [46]. Note, however, 
that our study of asymmetric salts near neutral surfaces 
reveals a closely related phenomenon: the generation 
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of an effective non-zero surface charge due to the 
unequal ionic response to a neutral dielectric interface 
for asymmetric salts. A further point that possesses 
experimental relevance is the role played by surface 
charge inhomogeneity. Strong-coupling calculations 
show that an inhomogeneous surface charge distribution 
characterized by a vanishing average value gives rise to 
an attraction of ions towards the pore walls, but this 
effect disappears at the mean- field level [J? . For a better 
understanding of the limitations of the proposed model, 
more detailed comparison with MC/MD simulations are 
in order [4"81 149] . Finally, dynamical hindered transport 
effects [571 EH] sucn as hydrodynamic forces deserve to be 
properly included in the theory for practical applications. 
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Appendix A: Variational free energy 

For planar geometries (charged planes), the transla- 
tional invariance parallel to the plane, allows us to sig- 
nificantly simplify the problem by introducing the partial 
Fourier-transformation of the trial potential in the form 

dk ik-(r 



WD&^ru-rfiH J —e^-^v^z 1 ,^). (Al) 



By injecting the Fourier decomposition ( |Al[ ) into 
Eq. (17), the DH equation becomes 

~lT e ( z )ir + e ( z )[ fc2 + K l( z )]\vo 0*,2',k; k v {z)) 
oz oz 



where the first term in the integral follows from Fq and 
the second term from (Hq) . Finally, the unscreened Van 
der Waals contribution, which comes from the unscreened 
part of Fq , is given by 
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The technical details of the computation of F 3 can be 
found in Ref. [5U] ■ The last term of Eq. ( A5 1 simply 



corresponds to the free energy of a bulk electrolyte with 
a dielectric constant e w . In the above relations, S stands 
for the lateral area of the system. The dummy "charg- 
ing" parameter £ is usually introduced to compute the 
Debye-Hiickel free energy [5T]. It multiplies the Debye 
lengths of vq (z, z, k; K v (z)) in Eq. (A4| and the dielec- 
tric permittivities contained in the thermal average of 
the gradient in Eq. ( A5 ) . This later is defined as 



my 



dk 



(A6) 



(k 2 + d z d z ,) v c [z,z',k;^(z)]| 2 



where we have introduced £j (z) = i^+S, \_&b {z) — ig 1 ] 
and v c [z, z' , k; £^(z)] stands for the Fourier transformed 
Coulomb operator given by Eq. ^ with Bjerrum length 



i^(z). The quantity F3 defined in Eq. (A5| does not de- 



pend on the inverse screening length k v . Moreover, in 
order to satisfy the elect roneutrality, <po(z) must be con- 
stant in the salt-free parts of the system where £b(z) ^ 
Ib- Hence, F 3 does not depend on the potential <j>o(z). 



Appendix B: Variational choice for the neutral 
dielectric interface 
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The translational symmetry of the system enables us to 
express any thermodynamic quantity in terms of the par- 
tially Fourier-transformed Green's function Vo{z, z, k). 
The average electrostatic potential contribution to F v 
that follows from the average (H} reads 



S dz 



8tt£ r 



+ Ps(z)Mz) 



(A3) 



the kernel part is 
S 



dz 



Ib(z) 



(A4) 



V (z,Z,\t-,K v [z)^ft} - V (z,Z,k- 1 K v (z)) 



We report in this appendix the restricted variational 
piecewise 4>q(z) for a neutral dielectric interface which is 
a solution of 



d 2 <f>o 
dz 2 

dz 2 



for z < a 



for z > a 



(Bl) 
(B2) 



where 4>q(z) in both regions is joined by the continuity 
conditions 



4>o(a) = 4>o (a), 



dcf>< 



dz 



dz 



(B3) 



We also tried to introduce different variational screening 
lengths in the second term of the lhs. and in the rhs. 
of Eq. (B2| without any significant improvement at the 



variational level. For this reason, we opted for a single 
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inverse variational screening length, k^. The solution of 
Eqs. pll)-p2| is 



ip for z < a, 

ip [1 + n^z — a)] e~ K *( z ~ a ) for z > a. 



where the coefficient c disappears when we impose the 



boundary and continuity conditions, Eq. (B3). The re- 



maining variational parameters are the constant poten- 
tial tp, the distance a and the inverse screening length 
K^. By substituting Eq. (B4) into Eq. (A3), we obtain 



the variational grand potential 



F V = V^- 
24tt 



8 
32^ 



A 2 K <p 2 



— Sp_ 



dz 



e -~w(z)+q-4> (z) _|_ 1- c -^w{z)-q + <t> {z) I 



whose solution is given by 



-^- + 2n(2-h) for z<h, 



Variational parameters introduced in this case are 
h, k^, and the charge renormalization 77, which 
takes into account non-linearities at the mean-field 
level jHj. The variational grand potential reads 

Fx = 2r ] (l + hk <p ) -r ] 2 {l/2 + hk 4> ) 
S 



2ir^k 



/ dze~ 9™ W cosh $ (z). (C6) 



In both cases, the boundary condition satisfied by 4>q 
is the Gauss law 



Appendix C: Variational choice for the charged 
dielectric interface 



dz 



2r] 



(C7) 



The two types of piecewise variational functions used 
for single charged surfaces are reported below. 

• The first trial potential obeys the salt-free equa- 
tion in the first zone and the NLPB solution in the 
second zone, 



where ?y = 1 for the non-linear case. It is important to 
stress that in the case of a charged interface, Eq. (C7) 



holds even if e ^ 0. In fact, since the left half-space is 
ion-free, 4>o(z) must be constant for z < in order to 
satisfy the global electroneutrality in the system. 



d 2 $ L 

dz 2 
d 2 $ h 



2S(z) for z<h, 



L sinh tfio = for z > h, 



dz 2 

whose solution is 

[ 4arctanh7 + 2(z — h) for z < h, 

lNL/;\ I 

90 \ z ) 



4arctanh(7e K ^ z for z > h, 



where 7 = k^ — w 1 



(CI) 



(C2) 



0. Variational parameters 

are h and k^, and the electrostatic contribution of 
the variational grand potential Eq. ( A3 1 is 



s 



h + 7 — 4arctanh7 k 2 
2tt S 47rS 



dze 



cosh^ L . 
(C3) 



Appendix D: Definition of the special functions 

The definition of the four special functions used in this 
work are reported below. 

Li »w = E£t> ew=Li„(i) (di) 



k>l 



P{x;y,z) 



f dtt v - 1 (l-t)'- 1 (D2) 
Jo 



2 F 1 (a,b;c;x) = }Xa)k{b)k(c)k 



k>0 



x 



(D3) 



where (a)/- = a!/(a — fc)!. 



Appendix E: Disjoining pressure for the neutral pore 



• The second type of trial potential obeys the salt- 
free equation with a charge renormalization in the 
first zone and the linearized Poisson-Boltzmann so- 
lution in the second zone, 



d 2 ^ 
dz 2 
d 2 4> L 
dz 2 



2rjS(z) for z < h, 
k\§\ = for z>h, 



(C4) 



The net pressure between plates is defined as 



1 8F V ( k\ 



(El) 



where the subtracted term on the rhs. is the pressure 
of the bulk electrolyte. The total van der Waals free 
energy, which is simply the zeroth order contribution F 
to the variational grand potential Eq. ([9]), is with the 
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FIG. 17: (color online) Difference between the pressure and 
the screened van der Waals contribution vs cL/Ib for KblB = 
0.5, 1 and 1.5, from left to right (e = 0). 



constraint n v — Kb (there is no renormalization of the 



inverse screening length at this order), 



vdW 



d4 

24?r 8ird 



Li 2 (e~ 



1 



16nd 2 



Li 3 (e 



-2dKb\ 



and 



vdW 



1 9F v dw 
S dd 



24tt 



(E2) 



(E3) 



We illustrate in Fig. [17] the difference between the van 
der Waals pressure and the prediction of the variational 
calculation for k^Ib = 0.5, 1 and 1.5. We notice that the 
prediction of our variational calculation yields a very sim- 
ilar behavior to that illustrated in Fig. 8 of Ref. [3T] . The 
origin of the extra-attraction that follows from the vari- 
ational calculation was discussed in detail in the same 
article. This effect originates from the important ionic 
exclusion between the plates at small interplate separa- 
tion, an effect that can be captured within the variational 
approach. 
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